rm(list=ls())
setwd("/Volumes/Disk2/Future_PM/Dataverse")
library(fields); library(maps); library(ncdf4);library(abind)
source('Function/get_geo.R')
source('Function/get_met.R')
source('Function/functions.R')

#########
mai=c(0.1, 0.05, 0.2, 0.05)
mgp=c(1.2, 0.4, 0)
tcl=-0.2
ps=12

ss=load("Data/ACCMIP/mask_high_resolution.Rdata")
lon.in=lon.out
lat.in=lat.out

plot.slope=function(slopes1,lon.out,lat.out,sig1){
mask=sp.dissolve(land.mask,lon.in,lat.in, lon.out, lat.out)
mask[is.na(mask)]=FALSE
slopes1[mask<0.3]=NA
slopes1[slopes1>=2.5]=2.5;slopes1[slopes1<=-2.5]=-2.5
slopes1[sig1>0.05]=NA
ind1=(lon.out>=-130 & lon.out<=-62.5)
ind2=(lat.out>=25 & lat.out<=51)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
image(lon.out[ind1],lat.out[ind2], slopes1[ind1,ind2],zlim=c(-2.5,2.5),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='');
	map("world", add = T)
	box()		
	spdata= slopes1
	spdata[is.na(spdata)]=0
	if(max(spdata,na.rm=TRUE)>=1){
		clines=contourLines(lon.out, lat.out, spdata,level=c(1))
		for(k in 1:length(clines)){
		lines(clines[[k]]$x,clines[[k]]$y,lwd=1,col=1,lty=2)
		}
	}		
}

dev.new(width=6,height=5)
close.screen(all.screens = TRUE)
m <- rbind(c(0,1,0.9,1),
		c(0,0.33,0.6,0.9),c(0.33,0.67,0.6,0.9),c(0.67, 1,0.6,0.9),
		c(0,0.33,0.3,0.6), c(0.33,0.67,0.3,0.6), c(0.67,1,0.3,0.6),			  		  
        c(0.2,0.8,0.23,0.53))
split.screen(m)
screen(1)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
text(x=0.5,y=0.2,bquote('JJA slopes of '*PM[2.5]~ 'and temperature'),cex=1.2)

screen(2)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
ss=load("Data/ACCMIP/AQS_species_slopes.Rdata")
slopes1=total_slope
lon.out=lon.frame;lat.out=lat.frame
plot.slope(slopes1,lon.out,lat.out,total_sig)
title("(a) obs",font.main=1,cex.main=1)

screen(3)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
ap=load("Data/ACCMIP/CAM3.5_JJA-PM25-T.Rdata")
lat.out=abind(lat.out, max(lat.out)+diff(lat.out)[1])
slopes1=cbind(slopes1,array(NA,length(lon.out)))
plot.slope(slopes1,lon.out,lat.out,sig1)
title("(b) CAM3.5",font.main=1, cex.main =1)

screen(4)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
ap=load("Data/ACCMIP/GFDL-AM3_JJA-PM25-T.Rdata")
plot.slope(slopes1,lon.out,lat.out,sig1)
title("(c) GFDL-AM3",font.main=1, cex.main =1)

screen(5)
ap=load("Data/ACCMIP/GISS_Lee.Rdata")
plot.slope(slopes1,lon.out,lat.out,sig1)
title("(d) GISS",font.main=1, cex.main =1)

screen(6)
ap=load("Data/ACCMIP/MIROC-CHEM_JJA-PM25-T.Rdata")
plot.slope(slopes1,lon.out,lat.out,sig1)
title("(e) MIROC-CHEM",font.main=1, cex.main =1)

screen(7)
ss=load("Data/GEOS-Chem/GEOS-Chem_species_slopes.Rdata")
plot.slope(total_slope,lon.out,lat.out,total_sig)
title("(f) GEOS-Chem",font.main=1, cex.main =1)

screen(8)
mai=c(0.1, 0.1, 0.05, 0.2)
par(mai=mai, mgp=mgp, tcl=tcl, ps=9)
image.plot(legend.shrink=0.9,legend.only=TRUE,zlim=c(-2.5,2.5),col=rwb.colors(32),horizontal=TRUE,legend.width=2.5,axis.args = list(mgp = c(0, 0.5, 0),tcl=-0.2,padj=-1
))
text(bquote(''~mu*'g' ~m^-3*K^-1*''),x=par("usr")[2]-0.23,y=par("usr")[3]+0.27,srt=0, adj = 0,xpd=TRUE,cex=1.3)
